One-dimensional cluster growth and branching gels in colloidal systems with 
short-range depletion attraction and screened electrostatic repulsion 
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We report extensive numerical simulations of a simple model for charged colloidal particles in 
suspension with small non-adsorbing polymers. The chosen effective one-component interaction 
potential is composed of a short-range attractive part complemented by a Yukawa repulsive tail. 
We focus on the case where the screening length is comparable to the particle radius. Under these 
conditions, at low temperature, particles locally cluster into quasi one-dimensional aggregates which, 
via a branching mechanism, form a macroscopic percolating gel structure. We discuss gel formation 
and contrast it with the case of longer screening lengths, for which previous studies have shown that 
arrest is driven by the approach to a Yukawa glass of spherical clusters. We compare our results with 
recent experimental work on charged colloidal suspensions [A. I. Campbell et al. cond-mat/0412108 
Phys. Rev. Lett, in press]. 



Recent years have witnessed a progressive interest in 
the role of the inter-particle potential on controlling 
structure and dynamics of colloidal dispersions. Experi- 
ments 0, H H HM J^J^JjEIpL theory [H Gj, 03 
and simulation |15l Ho . Il7l fl8L Il9t l20j | studies have pro- 
vided evidence that when the hard-core repulsion is com- 
plemented simultaneously by a short range attraction (of 
finite depth) and by a screened electrostatic repulsion, 
particles tend to form aggregates, whose shape and size is 
sensitively de pen dent on the balance between attraction 
and repulsion [la. l2ll 122 . 123. l24j . In some cases, the sys- 
tem shows an equilibrium cluster phase, where particles 
associate and dissociate reversibly into clusters 0, 0, @ . 
Interestingly enough, these cluster phases appear not 
only in colloidal systems but also in proteins solutions, 
in the limit of low salt concentration^ |g, Estimates 
of the ground state configuration of isolated clusters of 
different size^ suggest that, when the clusters diam- 
eter exceeds the screening length, the shape of the ag- 
gregates crosses from spherical to linear. Evidence has 
been reported that, for appropriate tuning of the exter- 
nal control parameters, colloidal cluster phas es p rogres- 
sively evolve toward an arrested state[Hl5lli^.ll0j]. It has 
been suggested, and supported by numerical simulations, 
that, in the case of relatively large screening length (i.e. 
the case of preferentially spherical clusters), dynamic ar- 
rest may proceed via a glass transition mechanism, where 
clusters, acting as super-particles interacting via a renor- 
malized Yukawa potential, become confined by the re- 
pulsions created by their neighboring clusters [15|. This 
mechanism is, in all respects, identical to the glass transi- 
tion of Yukawa particles |2i||2(||22],|28| and leads, favored 
by the intrinsic polydispersity of the clusters induced by 
the growth process, to the realization of a Wigner glass. 
The simulation study^f| showed that the resulting ar- 
rested state is not percolating, i.e. the arrest transition 
can not be interpreted in terms of the formation of a 
bonded network of particles. 
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very recent experimental work[T(il] has reported 
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FIG. 1: A pictorial view of the Bernal spiral. Particles have 
been differently colored to highlight the presence of three 
strands. In this geometry, each particle has exactly six nearest 
neighbors. 



idcncc of arrest via linear cluster growth followed by 
percolation, in a system of charged colloidal particles. 
In the studied system, the short-range attraction, in- 
duced via depletion mechanism, is complemented by an 
electrostatic repulsion, with a Debye screening length £ 
estimated of the order of £/<r rs 0.65, where a indi- 
cates the hard core diameter of the colloidal particle. 
The quasi one-dimensional clusters observed via confo- 
cal microscopy are locally characterized by a Bernal- 
spiral geometry 29], the same structure found as clus- 
ter ground state configuration for the case of screening 
lengths smaller than a |lj| ■ The Bernal spiral, shown in 
Fig. ^ is composed of face sharing tetrahedra, in which 
each particle is connected to six neighbors. 

In this work we numerically investigate the possibility 
that, when the potential parameters are such that the 
Bernal spiral is the ground state structure for isolated 
clusters, macroscopic gels can be formed at large, but fi- 
nite, attraction strength, via a mechanism of branching 
favored by the small but finite thermal contributions. We 
explore the low packing fraction region for several values 
of the attractive interaction strength, to highlight the 
collective effects arising from cluster-cluster interactions 
and to assess under which external conditions, ground 
state predictions are valid. We carry our study along two 
routes. In both cases, we study a colloid-polymer mixture 
in the effective one-component description, i.e. assuming 
that the polymer size is much smaller than that of the 
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colloids. In the first route, we control the attraction be- 
tween colloidal particles via a temperature scale. In the 
second route — designed to make direct contact with the 
experimental work reported in Ref.|lOl] — we study an 
isothermal system where the repulsive part of the poten- 
tial is fixed, while the attractive part of the potential 
is varied according to the concentration of depletant, to 
model the strength of the polymer-induced depletion in- 
teraction. We will refer to these two set of simulations 
respectively as temperature and polymer concentration 
route, naming them after the respective relevant control 
parameters. 

We study — as a function of the packing fraction and 
of the attraction strength — the shape of the clusters 
(quantified via their fractal dimension) , the local geome- 
try around each particle, the inter-particle structure fac- 
tor, the connectivity properties of the system. We com- 
plement the static picture with information on the inter- 
particle bond lifetime and on the dynamics of self and 
collective properties. We compare these quantities for 
the two routes, and show that the two approaches pro- 
vide a similar description of cluster growth, percolation 
and gel formation. 



I. SIMULATION DETAILS 

We study a system composed of N = 2500 colloidal 
particles of diameter a and mass to in a cubic box of size 
L, as a function of the packing fraction 4> c = irpa 3 /6, 
where p = N/L 3 is the number density, and of the tem- 
perature T. The particles interact simultaneously via a 
short-range potential Vsr and a screened electrostatic 
repulsive interaction Vy. The short-range attraction is 
modeled for simplicity with the generalization to a = 18 
of the Lennard-Jones 2a — a potential, as proposed by 
Vliegenthart et al. 



V SR (r) = 4e 



(1) 



where e is the depth of the potential. The parameters 
a and e are chosen as units of length and energy respec- 
tively. We also consider kg = 1. For this choice of a 
the width of the attraction range is roughly 0.2a. The 
phase diagram of Vsn(r) has been studied previously [3(j 
and it is characterized by a rather flat gas-liquid coex- 
istence line, with a critical point located at — 0.43 
and cj) S c R ~ 0.225. 

The repulsive interaction is modeled by a Yukawa po- 
tential 



V Y (r) = A 



77T' 



(2) 



characterized by an amplitude A and a screening length 
£. We focus on the case £ = 0.5c and A = 8e, for 
which the minimum of the pair potential Vsr + Vy is 



PV(r) 

r -10 




FIG. 2: Interaction potential 0V(r) = P[V SR (r) + Vy(r)] for 
different values of the polymer concentration <\> v . Here A = 
lOe, i = 0.65a, a = 10 and (3e = -U.0<j> p ^. 



located at r = 1.042c corresponding to a potential en- 
ergy E m i n = — 0.52e. With the present choice of A and 
£, the ground state configuration of an isolated cluster 
is known to be the one-dimensional Bernal spiral, shown 
in Fig. IT f Tg | . The Bernal spiral structure is composed 
of face-sharing tetrahedrons, resulting in three twisting 
strands of particles in such a way that each particle has 
six nearest neighbors. In this geometry, for large N, the 
potential energy per particle E is E — —1.36 + 2.10/iV 
(always in units of e). In the bulk of the spiral (far from 
side effects) E is about three times E m i n , confirming that 
the attractive interaction with the six neighbors provides 
most of the binding energy. 

In parallel, we also study the case in which the magni- 
tude of the attractive part changes to mimic the depen- 
dence of the depletion interaction on polymer concentra- 
tion <f>p. As in Ref. we choose e/ksT = — 14^ p , 
i.e. the attraction strength is assumed to depend lin- 
early on the fraction of free volume occupied by poly- 
mers <f>p. In this case, as in experiments, T is kept con- 
stant to UbT =1. To study a model as close as pos- 
sible to the experimental work of Ref. nj, we select 
£ = 0.65cr, A = 10e|3l|. and a = 10. For this value 
of a, r m i n = 1.07(7, in agreement with the position of 
the maximum in the radial distribution function <?(r), as 
extracted from data reported in Ref. [lfjj . 

The dependence of the potential shape with <f> p is 
shown in Fig[3 Note that V(r) changes from monoton- 
ically repulsive to repulsive with a local minimum (with 
V(r min ) > 0). Finally, for 4> p > 0.50, V(r) develops an 
attractive global minimum followed by a repulsive tail. 

In the rest of the present work, we will name T-route 
and (j)p route the two parallel set of simulations. The 
short-range nature of Vsr favors a very effective way to 
define pairs of bonded particles. Indeed, the resulting 
potential V(r) — Vsr + Vy has a well defined maximum 
located approximatively where the short range attrac- 
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FIG. 3: State points studied in this work, in the T — cf> c and 
4> P — 4>c planes, respectively for the T (circles) and <j> p (tri- 
angles) routes. Full symbols indicate state points where the 
equilibrium structure presents a spanning network of bonded 
particles. 



tion becomes negligible. In the following we will consider 
bonded (or nearest neighbors) all pairs of particles whose 
relative distance d < r max . In the T- route case, we fix 
Tmax = 1.28cr, while in the <p p case the bond distance be- 
tween two particles can be conveniently defined as lAa. 

In all simulations, time is measured in units of 
y/ ma 2 /e. For numerical reasons, the repulsive poten- 
tial is cut at r c = 8£, such that V Y (r c ) « 4.2 10~ 5 A 
All simulated state points are shown in Fig|3J Equilibra- 
tion is achieved with Newtonian dynamics, followed by 
a Brownian dynamics simulation, based on the scheme 
of Ref . |32j| . to produce equilibrium trajectories. In the 
case of Newtonian dynamics, the equation of motion have 
been integrated with a time step of At = 0.02. In the case 
of Brownian dynamics, At = 0.05 with a bare diffusion 
coefficient D Q = 0.005. Equilibration runs required, at 
the slowest states, more than 10 9 integration time steps, 
corresponding to about three months of computer time 
on a 1.6 GHz Pentium processor. 



II. EQUILIBRATION 

Simulations are started from high T (or correspond- 
ingly <j)p — 0) equilibrium configurations and quenched to 
the selected final state. During equilibration, a Berend- 
sen thermostat with a time constant of 10 is active, to 
dissipate the energy released in the clustering process. 
Following the quench, the time evolution of the potential 
energy E shows a significant drop. Equilibration becomes 
slower the lower the final T or the larger <j> p . It also slows 
down on lowering the colloid packing fraction <p c . The 
evolution of E following a quench is shown in Fig. 0] for 
the case <p c = 0.16. Around T < 0.07, equilibration can- 
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FIG. 4: Time dependence of the potential energy following 
a quench starting from high temperature (T = 1.0) for <j> c = 
0.16, for the T-route case. 



not be achieved within the simulation time and dynamic 
arrest takes place. In these conditions, extremely slow 
(logarithmic in time) drift of E is still present at long 
times. To provide evidence that equilibrium is reached 
during the Newtonian simulation we check that E is in- 
dependent on the previous history and that clusters re- 
versibly breaks and reform on changing T or <p c . Similar 
results are obtained following the </> p -path. 

The equilibration process is characterised by the pro- 
gressive formation of bonds between particles and the 
corresponding growth of particle's aggregates, named 
clusters. 

A quantification of the evolution of the structure of 
the system during equilibration can be provided by the 
structure factor S(q), defined as 



S(q) =< 



N ^ 



e -iq(ri-rj) > 



(3) 



where r\ indicates the coordinates of particle i. The 
S(q) evolution, shown in Fig. [S] is reminiscent of the 
initial stages of spinodal decomposition, showing a low q 
peak which grows in amplitude and moves to smaller and 
smaller q vectors. While in spinodal decomposition, the 
coarsening process proceeds endless, in the present case 
the evolution of the small q peak stops when equilibrium 
is reached. The presence of the low g-peak in S(q), at 
a finite wavevector, highlights the presence of an addi- 
tional characteristic length scale in the system, discussed 
in more details in the next section. 

Fig. HO shows the evolution of the shape of the largest 
cluster for the case (j> c = 0.125 and T — 0.08, one of the 
cases in which the average cluster size grows monotoni- 
cally in time. It is interesting to observe that at short 
times, the shape of the larger cluster is rather ramified, 
the potential energy is still large and locally the struc- 
ture is still very different from the six-coordinated ground 
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equilibration at cj> c — 0.125 and T = 0.08. 

FIG. 7: T (upper panel) and <f> p (lower panel) dependence 
of the normalized potential energy per particle E/E m i n at 
different <j>c values. E min — — 0.52e in the T-route case, 
while it depends on <f> p as E min — — 14(j>pkT + W(f m m), with 
r-min = 1-07 in the <f> p case. The corresponding value for the 
Bernal spiral configuration is also reported. 



T=0.08, <|> o=0.125 




pulsion enters into play, forcing thereby the system to re- 
arrange into the expected local configuration. This com- 
petition results also in a non-monotonic evolution, during 
the equilibration, of the mean cluster size, at some state 
points. 



t=200 t=600 t=20000 

III. EQUILIBRIUM PROPERTIES: STATICS 



FIG. 6: Snapshots of the largest cluster at three different 
times during the equilibration process. Here <j> c — 0.125 and 
T — 0.08. The cluster size is 72, 605 and 908, respectively at 
t = 200, t = 600 and t = 20000. 



state structure. Cluster arms are essentially composed 
by particles arranged along lines. At longer times, the 
cluster arms get thicker and thicker, and the local con- 
figuration approaches the characteristic one of the Bernal 
spiral, even if some parts of the original branching points 
persist in the final structure favoring the formation of a 
gel network. The evolution of the shape, complemented 
with the time dependence of E, suggests that at large at- 
traction strengths (low T or large <f> p ), the equilibration 
process can be conceptually separated into two parts: an 
initial relaxation which is closely reminiscent of the one 
which would take place if the potential was purely at- 
tractive, followed by a second rearrangement which sets 
in only after the coordination number has become signif- 
icant. At this point the competition of the long range re- 



A. Potential Energy 

The upper panel of Fig. shows the T dependence of 
E/E m i n at the studied values of <f> c . Around T w 0.2, 
E becomes negative, suggesting that the short-range 
attractive interaction becomes relevant. For lower T, 
0.1 < T < 0.2, E drops significantly, quickly reaching 
below T — 0.1 a value compatible with the ground state 
Bernal spiral configuration (also shown), once the vibra- 
tional components are properly accounted for. A similar 
behaviour is observed for the <p p dependendence, shown 
in the lower panel of Fig. [3 In the studied c -range, the 
4> c dependence of E is rather weak, especially for large 
attraction strengths. 



B. Cluster Size Distribution 

In this section we examine the cluster size distribution, 
as it evolves with (f> c and T. Standard algorithms are 
used to partition particles into clusters of size s and to 
evaluate the cluster size distribution n s and its moments. 
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The first moment of the cluster size distribution 



< s >-- 



N/N s 



(4) 



is connected to the inverse of the number of clusters N s , 
while the second moment < S2 > provides a representa- 
tive measure of the average cluster size 



< s 2 > 



(5) 



We also examine the connectivity properties of the 
equilibrium configurations. Configurations are consid- 
ered percolating when, accounting for periodic boundary 
conditions, an infinite cluster is present. The boundary 
between a percolating and a non-percolating state point 
has been defined by the probability of observing infinite 
clusters in 50% of the configurations. To provide an es- 
timate of the percolation locus we report in Fig. |3| the 
state points which are percolating. We note that, at this 
level, percolation is a geometric measure and it does not 
provide any information on the lifetime of the percolat- 
ing cluster. Indeed, at 4> c = 0.125, percolation is present 
both at high T, where we observe geometric percolation 
of clusters with bonds of very short life-time, and at low 
T where the particles are connected by energetic bonds 
of very long life time, as discussed below. The competi- 
tion between geometric and energetic percolation results 
in a intermediate temperature window where the system 
does not percolate, i.e. in a re-entrant percolation locus. 

The cluster size distribution n(s) is shown in Fig. [S] 
At T = 0.2 (where E f» and hence no significant 
bonding exists) upon increasing </> c , the distribution pro- 
gressively develops a power-law dependence with an ex- 
ponent t, consistent with the random percolation value 
r ss —2.2(33, Percolation is reached when 0.125 < 
cf) < 0.16. At slightly lower T, i.e. T = 0.15, the pic- 
ture remains qualitatively similar, except for a hint of 
non-monotonic behavior, around s w 10 — 20. On further 
lowering T, the number of clusters of size s < 10 drops 
significantly, to eventually disappear at T = 0.07. These 
results are observed at all studied densities. 

To frame the results presented above, we recall infor- 
mation previously obtained in the study of the ground 
state energy of isolated cluster of different size 191. For 
a cluster size s < 10, the addition of a monomer to an 
existing cluster lowers the energy per particle, since the 
gain associated to the formation of an additional attrac- 
tive short-range bond is not yet compensated by the in- 
creased number of repulsive interactions. However, when 
clusters have grown sufficiently, for s > 10 — 20, the en- 
ergy driving force for growing is reduced, since the en- 
ergy per particle does not significantly depend any longer 
on the cluster size^tj. Isolated clusters results carry on 
to the interacting clusters case since the relatively small 
screening length does not produce a significant cluster- 
cluster interaction. Indeed, the effective cluster-cluster 
potential will be characterized, to a first approximation, 
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Cluster size distribution n s at several T. In each 



panel, the full line represents the function n s 



-2.2 



by the same £0], which is short as compared to the 
distance between clusters. 

Around T w 0.2, E(T = 0.2) w and hence energy 
can not be the driving force for clustering. Thus, it is 
not a surprise that close to percolation n(s) ~ s~ T with 
t consistent with the random percolation value |33l l34|. 
The disappearance of clusters of size s < 10, which starts 
to be visible for T ^ 0.1, signals the progressive role 
of energy in controlling clustering. At the lowest inves- 
tigated T, energy has taken over and all clusters are 
formed by energetically convenient configurations. In 
this respect, we can think of the low T system as a fluid 
composed of super-aggregates, providing an effective re- 
normalization of the concept of "monomer" in the fluid. 
The small cluster-cluster interaction energy may favor 
a re-establishment of the random percolation geometries 
and characteristic exponents, as discussed in the follow- 
ing. 

Fig. shows the T and <fi c dependence of the sec- 
ond moment of the distribution, the average cluster size 

< S2 >, defined in Eq|31 for all non-percolating state 
points. Apart from <fi — 0.16, where configurations are 
percolating already before the physics of the short-range 
bonding sets in, percolation at small packing fractions 
is not reached at all temperatures we are able to equi- 
librate. At 4> = 0.125, a non-monotonic dependence of 

< S2 > (T) is observed, which we interpret as a crossover 
from the "random" percolation observed at high T to the 
bond-driven percolation, which becomes dominant at low 
T. The 4> c dependence of < S2 > is shown in the bottom 
panel. At all T, a monotonic growth is observed. 
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FIG. 9: Temperature and <f> c dependence of the second mo- 
ment of the cluster size distribution < si >, for the T-route 



C. Cluster shape 



FIG. 10: Typical largest cluster at 4> c = 0.08 for four 
different T values: from top left to bottom right T = 



T=0.15 




T=0.07 



A pictorial description of the shape of the larger clus- 
ter observed in a typical configuration at <p c — 0.08 and 
4> c — 0.125 for different T is shown in Figs.llOland 1111 In 
both cases, a progressive change of shape of the largest 
cluster is observed on cooling. A close look to the figures 
shows that on cooling particles become locally tetrahe- 
drally coordinated and that the loose high T bonding 
progressively crosses to a one dimensional arrangement 
of tetrahedrons. At the lowest T, the clusters are com- 
posed by large segments of Bernal spiral structures joined 
in branching points, the latter providing the mechanism 
for network formation. 

To quantify the cluster shape we study the the clus- 
ter size dependence of the cluster radius of gyration R g , 
denned as 



Ra 



1 N 



1/2 



(6) 



FIG. 11: Typical largest cluster at 
0.15,0.1,0.07. 



= 0.125 for T = 



where RcAf are the center of mass coordinates. For frac- 



tal aggregates, R g ~ s 1 '^, where df indicates the fractal 
dimension. The observed behavior of the clusters shape 
is very different at high and low T. Figll'2l shows R g 
vs s for two representative state points at T = 0.15, 
close to percolation. The typical shape of the cluster at 
these two state points is reported in Fig. 1101 and Fig. 1111 
When the cluster size is greater than 20 monomers, the 
fractal dimension is consistent with the random percola- 
tion value in three dimensions (df — 2.52 |33l l3^j ). This 
value confirms that at high T, as discussed previously, 
the energetic of the bonds is negligible as compared to 
cntropic effects and the cluster size grows on increasing 
<j) Cl mostly due to the increase in the average number of 
particles with a relative distance less than r max . At low 
T, an interesting phenomenon occurs, shown in Fig. 1131 
The very small clusters (s < 10) are rather compact 
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FIG. 12: Size dependence of the cluster gyration radius at 
T — 0.15 for two values of <f> c . The dashed line provides a 
reference slope for the random percolation df value. 



FIG. 13: Size dependence of the cluster gyration radius at 
T = 0.1 and several cj> c . Lines provide reference slopes for 
different df values. 



and df w 3 and indeed, in this size interval, the energy 
per particle in the cluster decreases on increasing cluster 
size[l!j. For clusters with intermediate size 10 < s < 100, 
df w 1.25, supporting the preferential one-dimensional 
nature of the elementary aggregation process, driven by 
the repulsive part of the potential. This df value is ob- 
served for all equilibrium cluster phases in which clusters 
of size 10 < s < 100 are dominant, with a small trend to- 
ward smaller values for smaller T and 4> c . This small de- 
value provides further evidence that in this size interval 
growth is essentially uniaxial, and that clusters of size 
100 or less are essentially composed by pieces of Bernal 
spirals joined by few branching points[ljJ (see Fig.llOland 
Fig. lll|) . For larger s values, a crossover toward df ~ 2.52 
is observed. This crossover suggests that for larger clus- 
ters the one-dimensional bundles have branched a signif- 
icant number of times, generating clusters whose geom- 
etry is again controlled by random percolation features. 
Pieces of Bernal spirals act as building blocks connected 
at branching points in a random fashion. 



D. Structure Factor 

As discussed in Sec. [HI the clustering process and the 
residual repulsive interactions between different clusters 
produce an additional low q peak in S (q) , located well be- 
low the location of the nearest neighbor peak (qa rj 27t) . 
Fig. [21 and Fig. ^| show respectively the T and 4> c de- 
pendence of S(q), in equilibrium. Data refer to both per- 
colating and non-percolating state points. We observe no 
dependence of the position (either with T or <f> c ) of the 
nearest-neighbour peak, consistent with the presence of a 
deep minimum in the interaction potential, which defines 
quite sharply the interparticle distance. The amplitude 
of the nearest-neighbour peak grows on decreasing T or 
increasing <fi c . The location of the cluster-cluster peak 



shows a weak <fi c dependence, almost absent at T — 0.2 
(and higher T), but which becomes more relevant at very 
low T. We note that, on isothermally increasing tp c , the 
location of the peak does not change even when percola- 
tion is crossed. On the other hand, the T-dependence is 
significant and the location of the peak moves to smaller 
q on decreasing T, suggesting the establishment of longer 
correlation lengths. The T and <p c trends are quite sim- 
ilar to those recently observed in concentrated protein 
solutions at low ionic strength^. In particular, in that 
paper, the independence of the cluster peak position on 
4> c was interpreted as evidence of a linear dependence 
of the equilibrium cluster size with <f> c . Indeed, if clus- 
ters are assumed to be rather monodisperse in size and 
if the inverse of the peak position is assumed to be a 
measure of the inter-cluster distance, the number clus- 
ter density has also to be independent on </> c |^- It is 
worth stressing that, in one of the first papers addressing 
the possibility of equilibrium cluster phases in colloidal 
systems |l2t l35j . the same relation between equilibrium 
cluster size and <j) c was presented, although its validity 
was limited to the case of clusters of size significantly 
larger than the one observed experimentally in Ref.jj], 
hinting to a wider validity of the relation suggested in 
Ref. 0| ■ Here we note that the independence of the S(q) 
cluster peak position with cf> c holds from very small <j) c up 
to values well beyond percolation, where an interpreta- 
tion in terms of finite clusters relative distance is clearly 
not valid. In the present study (of non spherical clusters), 
we can access both S(q) and the cluster size distribution. 
We note that, as shown in FigOO the cluster size distribu- 
tion is not peaked around a typical value. Actually, the 
cluster size is significantly non-monodisperse, expecially 
close to percolation. We also note that neither < s± > 
nor < S2 > (see Fig.JSJ scale linearly with (j> c , despite the 
constant position of the low q peak in S(q). 

It would be relevant to understand how the parame- 
ters A and £ entering the potential (see Eas 11120 . control 




FIG. 14: Wavevector q dependence of S(q) at three different 
T (T — 0.07,0.1,0.2) for the T-route case. For each T, data 
at three different <j> c are reported ( <j>c — 0.08, 0.125, 0.16). 




FIG. 15: Wavevector q dependence of S(q) at three different 
<j> c (4> c = 0.08, 0.125, 0.16) for the T-route case. For each <j> c , 
data at several different T are reported. 



the position of the cluster peak and its T and (f> c depen- 
dence. In the case of spherical clusters, it was possible to 
associate the peak position to the average distance be- 
tween clusters, since no percolation was observed. This 
explanation is not fully satisfactory for the present model, 
since, as can be seen in Fig.^|for the case of T = 0.2, the 
location of the peak is clearly the same both in the non- 
percolating state <fi = 0.125 and in the percolating state 
4> = 0.16. A better understanding of the quantities con- 
trolling the peak position is requested. A first attempt 
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FIG. 16: Average number of neighbors < n > as a function 
of T for different <f> c values. Note that, for all cj> c , all curves 
approach the < n >= 6 value characteristic of the geometry 
of the Bernal spiral. 

in this direction has been recently presented [l4j ]. 

E. Local Order 

A simple and useful indicator of local order is pro- 
vided by the average number of nearest neighbors < n > 
and by the associated distribution of nearest neighbors 
P(ri), which counts the fraction of particles surrounded 
by n neighbors within r max . As shown in Fig ll6l < n > 
grows upon progressively lowering T, approaching, in a 
non monotonous way, a coordination number of six. 

Fig.ll7lshows the T evolution of the distribution P(n). 
Again, a clear preference for local geometries with about 
six neighbors is displayed at low T a condition which 
is hardly observed in other materials in which particle- 
particle interaction is spherically symmetric. The value 
< n >— 6 is consistent with a local geometry of face- 
sharing tetrahedra. 

Another useful indicator of local order, which enables 
us to effectively quantify the local structure, is provided 
by the so-called local orientation order parameters q~i m (i) 
defined as, 

1 ^ 

q lm (i) = — Yi m (fy ) (7) 

where is the set of bonded neighbors of a particle 
i. The unit vector fy specifies the orientation of the 
bond between particles i and j. In a given coordinate 
frame, the orientation of the unit vector fy uniquely de- 
termines the polar and azimuthal angles 9ij and 4>ij ■ The 
Yi m {Qij i 4>ij ) = Yim ) ar e the corresponding spherical 
harmonics. Rotationally invariant local properties can be 
constructed by appropriate combinations of the q~i m (i). 
In particular, local order in crystalline solids, liquids and 
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FIG. 17: Distribution of the number of neighbors P(n) for 
several T at <j> c = 0.125 (for the T-route case). 



colloidal gels, has been quantified focusing on 
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m=—l 



1/2 



(8) 



and 



wi(i) = wi(i)/ 



E 



3/2 



(9) 



with 



wi{i) 



- E 

mi , ?d,2 , m3 , 
mi+m 2 +m3 = 



z / / 

mi mi rriz 



qi mi {i)qi m2 {i)<lim 3 {i)- 



The distributions of the qi and ibi parameters provide 
a sensitive measure of the local environment and bond 
organization. For example, dimers are characterized by 
qi = 1, ?«4 = 0.13 and wq = —0.09. A local tetrahedral 
order is characterized by large negative values of wq, up 
to the value —0.17 for the icosahedron j3(J. For the per- 
fect Bernal spiral of Fig^ the orientational order param- 
eters are determined as q± = 0.224, q 6 = 0.654, w 4 = 0.08 
and We = —0.148. Fig. 1181 shows the q^, q6, W4 and wq 
distributions and how they evolve with decreasing tem- 
perature for <j) — 0.125. We note that, upon cooling, the 
progressive presence of dimers and small clusters disap- 
pears and the distributions evolve toward a limiting form 
which appears to be specific of the Bernal spiral type of 
cluster. At low T, and in particular below T = 0.1, 
all distributions peak close to the characteristic values 
of the Bernal spiral. The local orientation order param- 
eters have been evaluated in the confocal experimental 
work of Ref. There, it was shown that the experi- 

mental data are consistent with the Bernal geometry. In 
the analysis of the experimental data, the position of the 
particles in the perfect spiral geometry was subjected to 
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FIG. 18: Temperature dependence of the rotational invariant 
distributions P(qi) (top) and P(wi) (bottom) for I = 4 and 
I — 6 at 4> = 0.125. Arrows indicate the ideal Bernal spiral 
values. In the ideal spiral, the local surrounding of all particles 
is identical and hence the rotational invariant distributions are 
delta functions. 



some random displacements, to account for thermal fluc- 
tuations, possible intrinsic errors in the localization of the 
particles and polydispersity in size (and/or charge) in the 
samples. After this procedure, the sharp peaks displayed 
in Figs. El an d El disappear and smooth distributions 
are obtained, which well compare with the experimental 
data. 

Fig. 1191 shows that, at low T, the distributions appear 
to be insensitive to <j) c , in agreement with observations 
in [lOj and supporting once more that the local struc- 
ture around the majority of the particles is similar to the 
ground state structure provided by the Bernal spiral. 
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FIG. 19: Packing fraction dependence of the rotational in- 
variant distributions P(qi) and P(u)i) for I = 4 and I = 6 
at T = 0.07. Note that, at this low T, no ^-dependence is 
present. 
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FIG. 20: Averaged mean square displacement < r > for the 
T-route (top) and the p -route (bottom)in log-log scale. In 
the top panel, <j> c = 0.16, while in the bottom one <f> c = 0.15. 



IV. DYNAMICS AND GEL FORMATION 

In this section, we present results for the particle dy- 
namics as a function of cj) c and T (in the T-route), or <fi p 
(in the <p p route). As for the equilibrium data shown in 
the previous section, dynamical quantities are evaluated 
from trajectories generated according to a Brownian dy- 
namics. The mean square displacement, < r 2 (t) >, aver- 
aged over all particles and several starting times is shown 
in Fig. |20| for one specific <j) c value both for the T and the 
4> p routes. 

Beyond the ballistic region (which extends up to < 
r 2 (t) >< 10~ 3 er 2 ), particles enter in a diffusive regime, 
composed of two different processes. A short transient 
where the bare self-diffusion coefficient D Q , set by the 
Brownian algorithm, dominates and long-time region 
when particles feel the interparticle bonding. At high T, 
in the latter regime, particles diffuse almost freely, with 
a diffusion coefficient not very different from the bare 
self-diffusion D a value. Upon cooling, < r 2 (i) > progres- 
sively develops a plateau, more evident for T < 0.1, which 
reaches the value rj 4 ■ 10 -2 . If we look at the p -data, 
we observe a very similar behaviour, with a very simi- 
lar plateau which develops for <j) p > 0.9. These results 
signal that particles become tightly caged, with a local- 
ization length not very different from the one observed 
in the case of dynamic arrest in glass forming systems, 
although in the present case caging is much less resolved. 
Increasing the attraction strength, the long time limit of 
< r 2 {t) > remains proportional to t, but with a smaller 
and smaller coefficient. 

A global view of the T and <j) c dependence of the slow 
dynamics is shown in Fig. 1211 where the long time limit 
of < r 2 (t) > /6t, i.e. the self diffusion coefficient D, 



is reported. While at high T the diffusion coefficient 
approaches the bare self-diffusion coefficient, on cooling, 
in the same T interval in which a substantial bonding 
takes place, D drops several order of magnitudes, signal- 
ing a significant slowing down of the dynamics and the 
approach to a dynamically arrest state. The same be- 
haviour is evident for the <p p route, where D approaches 
a very small value for 4> p > 1.1. 

It is interesting to note that for (f> c — 0.125 and 
<f> c = 0.16, the T dependence of D is compatible with 
a power law, with exponent 7rj « 2.2, not very differ- 
ent from the typical values of 70 predicted by MCT for 
simple liquids. The case of 4> c = 0.16 is particularly in- 
teresting, since at all T, the instantaneous configuration 
of the system is percolating, providing a clear example of 
the difference between percolation and dynamic arrest. 
Vanishing of D is observed only at very low T, well be- 
low percolation. It is tempting to state that, when the 
cluster-cluster interaction is weak as in the present case, 
dynamic arrest always requires the establishment of a 
percolating network of attractive bonds, though this is 
not a sufficient condition since the bond lifetime should 
be significantly long. When repulsive cluster-cluster in- 
teractions are not negligible, arrest at low <j> c can be gen- 
erated in the absence of percolation^] via a Yukawa 
glass mechanism. 

Another important quantity to characterize dynamic 
arrest (particularly relevant for attraction-driven slow- 
ing down .37]), is the bond correlation function <j>B (£), 
defined as 

Mt) = <5>y(*)Mo)>/PVo)]. (11) 

i<j 

Here riij(t) is 1 if two particles are bonded and oth- 
erwise, while Nb(0) = (J2i<j n ij(®)) ^ s tne number of 



FIG. 21: Temperature dependence of the normalized dif- 
fusion coefficient D/D , for different <f> c values. The short 
and long dashed lines represent power law fits with exponent 
7d = 2.15 and 70 = 2.37 and dynamic critical temperatures 
T d = 0.084 and T d = 0.091 respectively for (j> c = 0.125 and 
4> c = 0.16. The inset shows the corresponding quantity for 
the (f> p -route. 



bonds at t = 0. The average is taken over several dif- 
ferent starting times. 4>b counts which fraction of bonds 
found at time t = are still present after time t, in- 
dependently from any breaking-reforming intermediate 
process. 

Figure [221 shows the evolution of 4>b (t) with T and 
(pp. When dynamics slows down, the shape of 4>B(t) is 
preserved at all T or <fi p . The shape can be modeled 
with high accuracy with a stretched exponential function 
Aexp(— (t/r)P), with stretching exponent (3 « 0.73. 

An estimate of the average bond lifetime tb can be 
defined as tb = r//?r(l//3), where r and (3 are calculated 
via stretched exponential fits and T is the Euler Gamma 
function. 

FigEHshows tb versus T. Analogous considerations to 
those reported above in discussing the T-dependence of D 
apply. Indeed, tb{T) is consistent with a power law with 
exponent 7 T varying between 3.5 and 4.0, larger than the 
one found for D(T), but with consistent predictions for 
the diverging T. 

We notice that at T — 0.07, dynamics is extremely 
slow and bonds are almost unbroken in the time win- 
dow explored in the simulation. It would be interest- 
ing to find out if the T dependence of tb crosses to a 
different functional form at low T when all bonds are 
formed and if such cross-over bears some analogies to the 
cross-over from power-law to super Arrhenius observed in 
glass forming molecular systems. Unfortunately, as in the 
molecular glass cases, the time scale today available to 
simulation studies does not allow us to resolve this issue. 

To further compare the arrest observed in the present 
system and the slowing down of the dynamics observed 



FIG. 22: Bond correlation function <j> B (t) for C = 0.16 (T- 
route, top), and for <j> c = 0.15 (0 p -route, bottom). The 4>B(t) 
shape can be well fitted by a stretched exponential function 
with stretching exponent j3 = 0.73 (dashed line superimposed 
to the T — 0.12 curve). 




T 

FIG. 23: Temperature dependence of the bond lifetime tb at 
all studied densities. The short and long dashed lines repre- 
sent power law fits with exponent 7 T — 3.5 and Jt — 4.0 and 
dynamic critical temperatures T d = 0.084 and T d = 0.085 re- 
spectively for (j> c = 0.125 and (j> c = 0.16. The T = 0.07 point, 
not included in the fits, is shown here only as an indication, 
since equilibrium is not properly reached at this T. 

in other systems close to dynamic arrest, we calculate 
the collective intermediate scattering function F(q, t), de- 
fined as, 

F{q,t) =< ±-J2e~ im{t) ~ Pim > (12) 

where the average is calculated over different starting ini- 
tial times. Fig. 1241 shows the g-dependence of the F(q,t) 
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FIG. 24: Wavevector q dependence of the intermediate scat- 
tering function F(q,t) at T = 0.07,0.10,0.12 ( from top to 
bottom ) for 4> c = 0.16. The reported qa values are respec- 
tively 0.33, 0.78, 1.56, 2.34, 3.12, 4.68. 



at three different T values. The decay of the correla- 
tion functions does not show any appreciable intermedi- 
ate plateau for any q. The functional form of the de- 
cay is strongly dependent on q, crossing from an almost 
log(t) decay at small q to a less stretched decay at large 
q values. At the lowest T (T = 0.07), F(q,t) does not 
decay to zero any longer, confirming that a non-ergodic 
state has been reached. The non-ergodic behaviour man- 
ifests for very small values of qa, in the range of the low 
q peak in S(q), while ergodicity is restored at nearest 
neighbor length. Fig. [23 contrasts, at fixed q value, the 
T-dependence of the dynamics. The shape of F(q,t) is 
sufficiently different to conclude that time-temperature 
superposition does not hold for this observable. We also 
note that at very low <f> c {4>c = 0.04 or 0.08) all density 
correlation functions decay to zero, within the explored 
time window, suggesting that cluster diffusion allows for 
the decay of density fluctuation, even in the presence of 
a non-ergodic bond restructuring process. This suggests 
that, at low <fi c , density fluctuations are ergodic in the 
absence of percolation. 



V. CONCLUSIONS: 

In this work we have presented a detailed analysis of 
the structural and dynamic properties of a colloidal dis- 
persion in which the short range attraction is comple- 
mented by a screened electrostatic repulsion. We have 
studied one specific choice of the parameters controlling 
the repulsive potential. In particular, we have chosen 
a screening length comparable to the radius of the col- 
loidal particles. For this screening length, a study[l| 
of the ground state structure of isolated clusters showed 
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FIG. 25: Temperature dependence of the intermediate scat- 
tering function F(q,t) at cj) c = 0.16 and qa = 0.78. The 
reported T are 0.07, 0.1, 0.12, 0.15, 0.2, 0.25, 0.3, 0.4. 



that the preferential local structure is composed by a 
one-dimensional sequence of face-shared tetrahedra, gen- 
crating a local six-coordinated structure and a Bcrnal 
spiral shape. 

The collective behavior of the system is very much 
influenced by the competition between attraction and 
repulsion, which in the present model sets in when T 
becomes smaller than 0.2 (in units of the depth of the 
attractive part). The relative location of the particles, 
which for T > 0.2 is mostly controlled by translational 
entropy, for T < 0.2 depends more and more on energetic 
factors. Between T = 0.2 and T = 0.1, the number of 
bonded pairs increases significantly , and the local struc- 
ture evolves progressively toward the six-coordinated one 
characteristic of the Bernal spiral. At the lowest studied 
T,T = 0.07, the cluster shape becomes independent of (f) c 
and the ground state local configuration becomes dom- 
inant. The cluster size distributions at low T show a 
very clear suppression of clusters of size < 10, the size 
requested for the establishment of a bulk component in 
the spiral configuration. 

Although the majority of particles tends to preferen- 
tially sit in the 6-coordinated configuration, some parti- 
cles are located in defective regions of the spiral, which 
act as branching points and favor the formation of large 
ramified fractal clusters, whose elementary units are spi- 
rals of finite size. It is interesting to investigate if the 
small energetic cost of branching allows us to model the 
spiral segments as re-normalized monomers. In support 
of this possibility, we have detected a progressive in- 
crease of the cluster fractal dimension for cluster of size 
s > 100. We have also shown that, consistent with the 
ground state calculations, clusters of size s < 10 are al- 
most spherical, while clusters of size 10 < s < 100 are 
characterized by df w 1.25. 

The one-dimensional growth followed by a dynamic ar- 
rest phenomenon, observed in this work is reminiscent of 
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the aggregation process in several protein solution sys- 
tems |38l 13^. liol lil) . In this class of protein solutions, 
a variation in the external control parameters (temper- 
ature, ionic strength, p.H.) often trigger an aggrega- 
tion process of proteins into cylindric clusters which, by 
branching mechanisms, form a macroscopic gel, similarly 
to what takes place in the system here investigated. Re- 
sults reported in this work confirm that, as speculated 
in Ref . [13 , there is a range of small but finite tempera- 
tures in which branching of the one dimensional structure 
is preferred to cluster breaking and that such branching 
does indeed help establishing a connected three dimen- 
sional network. 

It is important to stress that the dynamic arrest mech- 
anism observed in this work is very different from the one 
observed numerically for the case of £ « 1.2cr|l5|. In that 
case, clusters grow mostly spherical and do not present 
branching points. The slowing down of the dynamics 
in the £ « 1.2a case arises from the residual repulsive 
cluster cluster interaction, resulting in the formation of 
a cluster phase or a repulsive cluster glass, analogous to 
the mechanisms suggested for Wigner glass systems. In- 
deed, in the arrested state, no percolation was detected. 
The arrested state generated via a Wigner glass tran- 
sition discussed in Ref.0] and the one generated via 
branching of one-dimensional clusters discussed in this 
work, although differing only by modest changes in the 
experimental conditions, are probably characterized by 
significantly different visco-elastic properties. Indeed we 
expect that the Wigner glass will be much weaker than 
the stiff percolating structure generated by a continuous 
sequence of particles tightly bounded to six neighbors. 

The system studied in this work is a good candidate 
for a thorough comparison with the slowing down charac- 
teristic of glass forming materials. The numeric "exact" 
equilibrium particle structure factor could be used as in- 
put in the mode coupling theory, along the lines theoret- 
ically suggested in Ref. ,13j to provide a full comparison 
of the theoretical predictions for the arrest line as well 
as for the shape of the correlation functions. It would be 
interesting to quantify the role of the cluster pre-peak in 



the structure factor in the predicted slowing down of the 
dynamics. 

Results presented in this work also provide further ex- 
ample of the existence of equilibrium cluster phases, a 
phenomenon which is recently receiving a considerable 
interest. Cluster phases have recently been investigated 
in systems as different as pro tein solutions B, 0, |t| , col- 
loidal dispersions [H [3. [EL |lOl|. laponite l42j. liposomic 
solutions p. I43L [5 , l45l l46l] 47 | . star-polvmers[ll| , aqueous 
solutions of silver iodide |3j , metal oxides [jj and in recent 
numerical studies [T5L Il6t Il7i l2*cj . In all cases, the combi- 
nation of the repulsive interactions with the short range 
attraction appears to be crucial in stabilizing the cluster 
phase. The high sensitivity of the cluster shape and the 
final topology of the arrested state on the detailed bal- 
ance between range and amplitude of the attractive and 
repulsive part of the potential brought forward by this 
and previous studies add new challenges to the modern 
research in soft condensed matter and to the possible 
technological exploitations of these new materials. 

A final remark concerns the use of an effective poten- 
tial, with state-independent parameters for the descrip- 
tion of systems in which the screening length can be a 
function of the colloid packing fraction and in which the 
significant changes in structure with T (or concentration 
of depletant) may lead to relevant changes in the clus- 
ter surface potential or in the spatial distribution of ions. 
The similarity between the numerical data reported in 
this manuscript and the closely related experimental re- 
sults suggest that, despite the approximation adopted in 
the numerical work, the essence of the arrest phenomenon 
is captured by the present models. 
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